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A LEAST SQUARES METHOD FOR THE REDUCTION OF 


FREE -OSCILLATION DATA 

By Phillip R. Wilcox and William L. Crawford 
Ames Research Center 


SUMMARY 


The classical least squares curve fitting method is used to determine the 
frequency, amplitude, damping ratio, phase angle, and zero offset of both a 
one- and two-degree-of -freedom system from free oscillation data. 

The method is applied to a number of experimental transients with good 
results. Where possible, comparisons are made with the results of other 
methods. The least-squares method is found to be particularly useful in the 
analysis of two-degree-of -freedom systems where other techniques are difficult 
or impossible to apply. 


INTRODUCTION 


In wind-tunnel testing when the dynamic stability of a configuration is 
one of the parameters to be determined, either free -oscillation or forced- 
oscillation techniques can be used. The free -oscillation technique is easier 
to implement and was used recently at Ames Research Center in testing axi- 
symmetric hammerhead models (ref. l) . A typical model was mounted in the 
wind tunnel on a free -oscillation balance and a disturbance in the model atti- 
tude was introduced by rotating the model and then quickly releasing it. A 
continuous signal proportional to the model attitude was obtained and later 
passed through an analog to digital converter to provide discrete measurements 
at equal time intervals. When this record corresponds to a one-degree -of - 
freedom system, there are several methods commonly used to obtain the damping 
ratio and natural frequency. One of these methods (called here amplitude 
response method) is to use the amplitude response curve (ref. 2) as indicated 
in figure 1. This curve is obtained from the Fourier transform of the free- 
oscillation transient. Another method (called here peak amplitude method) is 
to plot the values of the peak amplitude on a semilog plot and then calculate 
the damping ratio as indicated in figure 2. 

During the reduction of the data from these tests it was found the 
methods mentioned above failed to yield consistent data for a number of tran- 
sients. These transients all showed evidence of multiple-mode oscillations. 

In an effort to find a data analysis technique that could be applied to tran- 
sients of this type, a least-squares method was developed. It is the purpose 
of this paper to describe this method which will analyze either a one- or two- 
degree-of -freedom system rapidly and accurately, and to describe how the 



results obtained by this method compare with the results obtained by the 
amplitude response and peak amplitude methods. 

NOTATION 

A amplitude of the envelope of curve 1 at t = 0, same units as y(t) 

B amplitude of the envelope of curve 2 at t = 0, same units as y(t) 

C zero offset, same units as y(t) 

f phase angle of curve 1, radians 

g phase angle of curve 2, radians 

n number of data points 

p natural frequency of curve 1, radians/sec 

q natural frequency of curve 2, radians/sec 

t time, sec 

y(t) amplitude of transient 

y (t) measured values of transient 

y c (t) calculated values of transient 
a -£iP, 1/ sec 

0 -izn* l/sec 

£ ! damping ratio of curve 1 

damping ratio of curve 2 

DESCRIPTION OF THE LEAST -SQUARES METHOD 

The response of the one- and two-degree-of -freedom transient is assumed 
to be of the form 

y(t ) = Ae at sin(pt + f) + C (l) 

and 

y(t ) = Ae at sin(pt + f) + Be^ sin(qt + g) + C (2) 
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The problem is to determine the constants in these equations so that the 
equation best describes, in a least -squares sense, a set of experimental data* 

Appendix A is a detailed description of the derivation of the equations 
needed to solve for the unknowns in equation (l) or (2) by the LSM (Least 
Squares Method) . The analysis results in a matrix equation that can be solved 
by normal matrix techniques* To apply this method to a two -degree -of -freedom 
system, first assign a vector (A 0 ,aojPo/£o>C 0 jBo,Po><lo jgo) of initial esti- 
mates for all the unknowns. Using these values, calculate both y (denoted by 
y c ) and the partial derivatives of y with respect to each of the unknowns 
(§y/dA,dy/d<x, . . * , dy/dg) at times corresponding to each of the n data 
points. After doing the indicated summations the following matrix equation 
for corrections to the initial estimates is obtained. 
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After solving these nine equations for the corrections (AA,Ax, . . . ,Ag) , 
they are added to the original estimates of the constants and the process is 
repeated. This iteration process is continued until the corrections are arbi- 
trarily small (assuming convergence) at which time the estimate vector contains 
values of the constants which best describe the experimental data. 

This procedure does not guarantee convergence. It has been found, how- 
ever, that if reasonable initial estimates can be made of the constants (in 
particular the p and q) , then convergence is assured and is very rapid. 


PRESENTATION OF RESULTS 


One Degree of Freedom 

The transients of the three models selected are plotted in figure 3 * 

They were chosen because they are representative of the range of damping 
ratios encountered during the actual testing program described in reference 1. 
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The peak amplitude and amplitude response plots of each transient ar- 
shown in figures 4 and 5 , respectively. The Share Library Program, AAHAfb, 
mode 3 (ref. 3 ) was used to calculate the amplitude response of the transier. 
(i.e., Fourier transforms). This program was converted to MAP for use with 
Fortran IV programs by the staff of the Computation and Analysis Branch at 
Ames. The damping ratios of the peak amplitude plots were calculated as indi- 
cated in figure 2 while the frequencies were determined from figure 3 by 
counting the number of cycles in a given time increment. The damping ratios 
and frequencies of the amplitude response plots were calculated as indicated 
in figure 1. These results, the results obtained from the LSM, and the nor- 
malized standard deviation, or SD, (see section on Error) are recorded in 
table I. From this table it can be seen that the SD’s, calculated with the 
parameters from the LSM, for the three transients are very small. Since 
these values are less than the accuracy of the data, the values obtained from 
the LSM for the frequencies and damping ratios are assumed to be correct. 

Now, with the parameters of the LSM as the basis, the errors of the other 
methods are as follows: 


Transient A 


Transient B 


Transient C 


Method 


£ L £ L 

- 0.36 $ +7.14 $ -2.09% -5. 

- 0 . 0 6% + 185 . 60 $ + 0 . 17 $ + 37 . 


£ L 

- 1 . 68 $ +2.17$ Peak amplitude 

+0.20$ +34.35$ Amplitude response 


Prom this, it can be seen that all three methods compare veil on frequency 
but the amplitude response method gives damping ratios that do not agree with 
those obtained from the other tvo methods. 

The solutions obtained from the LSM have been plotted against the experi- 
mental data in figure 6 . It can be seen from this plot that there is very 
good agreement betveen the generated curves and the transient data. 


Tvo Degree of Freedom 

The three transients chosen for the tvo -degree -of -freedom system have 
been plotted in figure 7* They vere chosen to give a variety of tvo -degree - 
of -freedom cases in vhich the applicability of peak amplitude and amplitude 
response methods becomes doubtful. The peak amplitude plots are as shovn 
in figure 8 , and the fact that no straight line can adequately represent the 
peak values is apparent. The amplitude response plots are shovn in figure 9* 
and although the frequencies are veil defined, the bandvidth measurement for 
the second mode of transient E vas uncertain. These results, along vith those 
obtained from the LSM, and the normalized standard deviation, or SD, are 
recorded in table II. Here again the parameters as calculated by the LSM are 
assumed to be correct for the reason stated previously. With these parameters 
as a basis, the errors for the amplitude response parameters are as follovs: 
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Transient F 



Transient D 

Transient E 

p 

+0.02$ 

+0.18$ 

a 

+0.53$ 

+3.68$ 


+ 18 . 45 $ 

-15.15$ 

£2 

-25.30$ 

— 


+ 0 . 33 $ 


+ 0 . 10 $ 


+33.80$ 

+22 • 10$ 


From this, it can he seen that although hoth methods agree well on frequency, 
they disagree on the damping ratios. 

The results of the LSM for transient D are shown in figure 10 where the 
individual transients, curves 1 and 2, and their sum have heen plotted. The 
sum is plotted against the experimental data. The results obtained for tran- 
sients E and F are plotted in figures 11 and 12, respectively. It can be 
seen that there is again good agreement between the experimental data and the 
theoretical curve for each of the transients. 


Error 

The present LSM has been programmed to assume that it has successfully 
arrived at the values of the parameters when the incremental change in all 
parameters is less than 0.05 percent of the previous value of that parameter. 

As an aid in judging how well the theoretical curve fits the experi- 
mental data, the normalized standard deviation is computed. 



X 100 


This parameter has been recorded in table I for the one -degree-of -freedom 
transients and in table II for the two -degree -of -freedom transients. The 
poorest fit was for transient D, where SD = 1.98. Thus, in all cases, very 
satisfactory fits were obtained. 

Since the second frequency for transient E was not pronounced (see 
fig. 9)> this transient was analyzed also as a one -degree -of -freedom system. 
The best fit obtainable under this assumption was SD = 2.11, significantly 
poorer than the two -degree -of -freedom case where SD = 1.29. 
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CONCLUSIONS 


A least squares method has been developed for the analysis of free- 
oseillation data and a comparison of this technique -with the peak amplitude 
and amplitude response methods has led to the folio-wing conclusions: 

1. For the one -degree -of -freedom transient data, all three methods give 
good answers for the frequency but the amplitude response method gives poor 
values for the damping ratio, 

2. For the two -degree -of -freedom transient data, the least squares 
method gave a good fit for cases where the other methods were not practical 
to apply. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif. 94035, Jan. 8 , 1968 
124-11-04-09-00-21 
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APPENDIX A 


MATHEMATICAL DERIVATION OF NUMERICAL SOLUTION 


It is assumed that the free -oscillation transients considered in this 
paper can he modeled as the sum of one or more damped sine waves; that is 


y(t) = Ae 00 ^ sin(pt + f) + Be^(qt + g) + . . . + C (Al) 

where t is the time variable; C is a zero offset correction; A, B, . . . 
are amplitude factors; a, (3, • • . are functions of the damping ratios; 
p, q, . • . are frequencies; and f, g, . • • are phase angles. 

The problem is to determine these parameters from the transient data 
such that equation (Al) best describes the transient data in the least squares 
sense. 


The usual method of using the least squares criterion (ref. 4) directly 
for fitting equation (Al) to the transient, y(t), would result in a set of 
nonlinear equations for which there is no known analytic solution. However, 
an iterative method of determining successive corrections to estimates of the 
parameters can be applied. This method (refs* 4 and 5) which results in 
linear equations is derived in the following paragraphs. 

The derivation will be made for the solution of the equation 

y(t) = Ae a ^ sin(pt + f) + Be^ sin(qt + g) + C (A2) 


Extending the results to any number of damped sine waves is straightforward. 
Equation (A2) can be written in the symbolic form 


y = y(t) = y(t,A,a,p,f,C,B,0,q,g) 


(A3) 


To apply the least squares method, we need to express the relationship 
between the empirical data and equation (A3) as a set of residual equations. 
Given n data points y^_, i = 1, 2, . . . , n, if Vj_ represents the differ- 
ence or residual between the data y-j_, taken at time tj_, and the value cal- 
culated from equation (A2) at time ti, we obtain the following residual 
equations : 
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Vi = y(ti,A,a, . . • , g) - yi 


a — y(t^,A,a, • • • > g) ~ y ± 


V n = y( t n > A > a > • 


• > g) • y r 


If we now assume that we can obtain approximate values A 0 ,ao,Po,Bo.>Co, 
Bo,3o><lo,go Bor the true parameters A, a, . . . , g, then we can express the 
true parameters as 


A = 

A 

+ 

AA 

B = 

B 

+ 

AB 


o 




o 



a =■ 

a o 

+ 

Ax 

3 = 

Po 

+ 

Ap 

P = 

P o 

+ 

Ap 

q. = 

% 

+ 

Aq 

f = 

f 

o 

+ 

Af 

g = 

g o 

+ 

Ag 

C = 

C 

+ 

AC 






where AA,Ax, . • . , Ag are corrections to be determined. Substituting 
equation (A5) in (A4) we obtain the residual equations 


= y(ti,A 0 + AA,a 0 + Ax, . . . , g Q + Ag) - yi 


r , = y(t, ,A + AA,a + Ax, • . . , g + Ag) - y. 


v = y(t -A 4- AA y a + Ax* # 
n J v n' o * o 7 


. , g Q + Ag) - y n 


For simplicity let us consider only the ith residual equation, which 
can he written as 
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(AT) 


V ± + y ± = y(t ± ,A o + AA,a o + Ax, . . . , g Q + Ag) 

By considering the right side of equation (A 7 ) as a function of A, a, . % . , g 
with t constant, we can expand it by Taylor 1 s theorem for a function of sev- 
eral variables (refs, 5 and 6) around the point (A 0 ,ooj • • • , go)* This 
expansion is valid since equation (A2) has continuous partial derivatives of 
all orders with respect to A, a, . . . , and g. 

Taylor’s expansion yields 


v. 

x 


y i 


y(t 


- jA ,a y 

X 7 o 7 o 7 



+ 



0=1 





3 

) y(t,A,a, 


9 


g) 




(A8) 


where * means that after differentiation the numerical values of the partials 
are calculated for t = ti, A - A 0 , a = ao; P=Po> f = foj C = C 0 j B = B 0 , 

3 ^ Poj q. = <lo> g = go- 

Expanding equation (a8 ) gives 

Vi + y ± = y(t ± ,A o ,a o , . . . , g Q ) 

oo 

+ M (Ml * ^ (Ml + • • • - ^ (Ml i r j (A9) 

0=2 

Since sin(pt + f, qt + g) ^ 1, cos(pt + f, qt + g) ^ 1 for all (pt + f, 
qt + g) ; and at, pt, A, and B are bounded; it is clear that the Rj consist 
of sums of terms Rjk the or ^ er 


R - order 
jk 


max(l,A,B)max(AA 2 ,AA Ax, 


AA Ag, 


Ax 2 , Ax Ap, . . • , Ap 2 , . • • , Af , • • » , Ag ) 


(A10) 


Since all the Aa, Ax, . . . , Ag are considerably less than 1, and all the 
partials are bounded, all the Rj approach zero as j gets large. This 
result, together with the fact that all the partials exist and are continuous 
in the range 0 £ t ^ 1 is a sufficient condition (ref. 7) for equations (A8) 
and (A9) to be valid. 
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If we now approximate y(ti,A 0 + AA,cxo + Ax, • 
the second and higher order terms of equation (A9) , 


go + Ag) hy dropping 


v. + y . = y(t. ,A ,<x , 
x 1 J i 7 o’ o’ 


g Q ) + AA 



+ Ax 



+ 



(All) 


or 


v. = AA 

i 



-f 



+ y(t.,A o , v 



(A12) 


If we let 7i = y(ti,A 0 ,ao, 


go) “ yi, equation (A12) becomes 



. . + Ag 



i = 1,2, 


n 

(A13) 


Equation (A13) now represents a set of n residual equations which are linear 
in the corrections AA, Ax, . . . , Ag and can thus be solved for the correc- 
tions by a least squares technique. 

The least squares criterion is that 

i=i 

This criterion will be satisfied when the first partials with respect to all 
the unknowns, evaluated at the point (A 0 ,a 0 , . . . , go) , are equal to zero. 
Thus for the partial with respect to AA: 


n 

v_. 2 will be a minimum (ref. 7)* 



i=i i=x 


Substituting equation (A13) for Vj_ becomes 


n 



i=i 




4 * • 


+ Ag 
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or 


n 



Distributing the summation and taking the constants A A, Ax> . . . , Ag 
outside the summations gives 



i=i 


(A14) 


Proceeding similarly for each of the other corrections Ax, Ap, . • • , Ag 
results in nine linear equations in the nine unknown corrections. Numerical 
values for the corrections can now be obtained by matrix techniques (ref. 8). 
If we arrange into an array the coefficients of the unknowns in the set of 
nine equations derived in the same manner as equation (Al^) we obtain the 
square matrix U of order 9 - 




I (A15) 

n , v 2 

VAA 

i=l 
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We now define a column vector V as the left hand members of the nine 
equations 


or 



(A1 6 ) 


We now wish to determine the vector X such that UX = V. Jordon's 
method (ref. 9) is used to reduce U to I , the identity matrix, through 
a series of elementary transformations, which when applied to V result 
in X. The vector X now represents the corrections: 


X = 




(A17) 
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Since this result was obtained by a truncation of the Taylor series 
expansion in equation (All) , X contains only approximate values for the 
corrections. We, therefore, form a new set of approximate values Ax, ai, 
• • • 9 gi by adding the computed corrections to the initial estimates 


At = A + AA 
- L o 

cci = a + Ax 
o 


(A18) 


gl = g Q + Ag J 


For notational convenience we now let Ax be called A 0 , cci be called 
ao; • • • , and g x be called g Q , and return to equation (A12). We now 
proceed as before to obtain an improved approximation to the true parameters 
A, a, • . • , g. This process of iteration is repeated until there is no 
change (to the degree of accuracy obtainable from the input data) in the 
parameters A, a, . . . , g in equation (Al8) in two subsequent iterations. 

By using the computer program based on this analysis, we have shown that 
this iterative process does converge rapidly under a specific range of condi- 
tions. The primary condition is that the initial estimate of the frequencies 
be close to the true values. The choice of first and last data points pre- 
sented to the program is not critical because of the parameters included for 
phase angles. However, including data containing extraneous forces or vibra- 
tions, or data with nonstationary parameters within the time period used will 
definitely limit or negate convergence to meaningful values. 




.. n 
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TABLE I.- COMPARISON OF RESULTS FOR ONE -DEGREE -OF -FREEDOM TRANSIENTS 


Transient 


Method 


Frequency . , . . . , , Phase angle 

f/ n \ Damping ratio Amplitude 

l p /at) ( a ) 


radians 


Normalized 

Zero offset standard deviation 
(c) (SD) 


Peak amplitude 

16.6 

Amplitude response 

16.65 

Least squares 

16.66 

Peak amplitude 

17.3 

Amplitude response 1 

17.70 

Least squares 

17.67 

Peak amplitude 

19.8 

Amplitude response 

20.30 

Least squares 

20.26 


- 0.0044 

0.20 

.0138 

.78 

.0257 

1.42 






TABLE II.- COMPARISON OF RESULTS FOR TWO -DEGREE -OF -FREEDOM TRANSIENTS 




Frequency 

Damping ratio 

Amplitude 

Phase angle 

Zero 

Normalized 

standard 

Transient 

Method 

Curve 1 

(p/2rt) 

Hz 

Curve 2 
(q./2it) 
Hz 

Curve 1 

(£i=-cc/p) 

Curve 2 
(^2 = "P/ 9.) 

Curve 1 
(A) 

Curve 2 
(B) 

Curve 1 

(f) 

radians 

Curve 2 
(g) 

radians 

offset 

(c) 

deviation 

(SB) 

D 

j Peak 
; amplitude. 


— 

— 

— 








: Amplitude 
' response 

25.50 

! 32.45 

i 

0.0088 

0.0157 








Least 

squares 

25.46 

1 

32.28 

.007 43 

.02097 

1.0161 

; 1.3805 

- 0.8603 

0.3161 

-0.0147 

1.98 

E 

Peak 

amplitude 

— 

— 

— 

— 








Amplitude 

response 

21.70 

. 23.55 

• 0235 

— 







! 

Least 

squares 

21.66 

24.45 

.02769 

.01546 

- 1.9502 

-.4288 

.0017 

2.499 

.0009 

1.29 

F 

Peak 

amplitude 

— 

— 

— 

— 








Amplitude 

response 

45.20 

102.20 

.0149 

.0069 








Least 

squares 

45.05 

102.10 

.01114 

. .00565 

.6614 

-1.7834 

.4784 

.0429 

-.0141 

1-57 
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Figure 6.- Results of the transients as calculated from the LSM 
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Figure 8.- Peak amplitude of transients used in the two-degree-of-freedom 

analysis . 
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(a) Curve I 



(b) Sum (curve I and curve 2) 



Figure 10,- Results of transient D as calculated from the LSM. 







National Aeronautics and Space Administration 
Washington, D. C. 20546 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS AND 
SPACE ADMINISTRATION 


OFFICIAL BUSINESS 


FIRST CLASS MAIL 


! 



i 


1 


1 1 / 


MASTER: 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 


"The aeronautical and space activities of the United States shall be • 
conducted so as to contribute ... to the expansion of human knotvl- 
edge ■ of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and Notes, 
and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



